This "flow downslope" example involves four sub-directories, layer, rho, sigma and z, in which the model is running in one of four coordinate configurations. To use this notebook it is assumed you have run each of those experiments in place and have kept the output generated in default configuration.

CAVEAT: This is a tutorial of how to make vertical section plots which also illustrates some poor ways to plot data for comparison. Read through to the end.

We will use matplotlib. The line %pylab inline loads all the necessary packages, including numpy and causes images to appear in the page. We will use scipy's netcdf package to read MOM6 output. Note that this only works if MOM6 is compiled with NETCDF=3.

To see this notebook with figures see https://gist.github.com/adcroft/dde8d3fafd77d0caaa5613e64f1d7eff.

Accessing MOM6 output

Let's open some MOM6 output files and see what variables we diagnosed. We'll use scipy.io.netcdf_file but you could also use the netCDF4 package and netCDF4.Dataset with the same syntax. Both functions return objects with members (functions or dictionaries) such as .variables:

The above output is the result of this diag_table:

We actually only asked for "e", "h", "u", "v", "salt" and "temp" in the diag_table. The other variables (such as "xh", "yq", etc) are all 1-dimensional "axis" variables associated with netcdf dimensions. This is a consequence of the CF convention used in MOM6 output.

Coordinates to use for plotting

Because MOM6 uses curvilinear coordinates in the horizontal, the axis variables for the horizontal dimensions are not always useful.

This particular experiment is in Cartesian coordinates so we can use the variables "xh" for the h-point locations and "xq" for the u-points. In this instance we can read the 1D axis variables to use as coordinates:

In addition to output requested in the diag_table, MOM6 normally writes a files with static coordinates and metrics to ocean_geometry.nc:

Here, you'll see some 2D coordinate variables such as geolon and geolonb that facilitates plotting when the coordinates are curvilinear. We can check that the data is Cartesian and matches the CF dimensions-variable we read earlier with

The above plot just confirms that we have the same values in numerous forms of the horizontal coordinates (for this experiment). In general, to make plan view plots for a 3D configuration you will always what to use the 2D coordinate data such as that in ocean_geometry.nc but here we will be making vertical section plots which only require the 1D form of horizontal coordinates.

Reading bottom depth from ocean_geometry.nc

Now let's plot the topography, which is contained in the variable D of ocean_geometry.nc.

"Depth" is a positive quantity for the ocean; the bottom is at height $z=-D(x,y)$. The topography visualized for this experiment is a flat shallow shelf on the left, a flat deep ocean on the right, and a linear slope in between.

The "flow downslope" example is ostensibly a 2D configuration in x-z. One would expect the y-dimension to be equal to 1 but it is instead equal to 4. Each of the 4 j-slices of an array will contain the exact same data. This is because a dimension in MOM6 can not be reduced to below the width of a parallelization halo width, which is typically equal to 4 or more.

This quickly illustrates that the model state is identical along the j-axis:

So from here on, we will use j=0 in all plots.

Exploring vertically distributed model output

Now let's look at some model data in multiple coordinate modes. We opened the output files above. The python variable layer_file is a handle to the netcdf file layer/prog.nc which is output using the traditional isopycnal (or stacked shallow water) mode of MOM6. The other models in this experiment are all ALE-mode emulating the z*-coordinate (python variable z_file for z/prog.nc), terrain-following sigma-coordinate (sigma_file for sigma/prog.nc), and continuous isopycnal coordinate (rho_file for rho/prog.nc).

The diagnosed variables in each of these modes was the same. However, some axis data changes meaning. For example, the vertical coordinate in layer mode is a "target density":

When the model is in ALE mode emulating a z* coordinate, then the vertical coordinate is height (although we report notional depth to aid ferret with plotting):

Let's look at salinity in the first record written by the model in each of these four coordinates. We'll plot the raw data without coordinates, i.e. in index space.

There is no topography apparent in the plots and the salinity structure is hard to make sense of! In layer mode there is not even any horizontal structure. This is because in layer mode density is homogeneous along a layer whereas in ALE mode density (in this case salinity) is allowed to vary along layers.

Plotting output from any model in index-space ignores the coordinates which determines where the data is physically located. The apparent absence of topography is a symptom of this. In MOM6, layers always contain data but layers have variable thickness which can even vanish. This is what thickness looks like for the above salinity panels:

The simplest distribution to explain is the terrain-following $\sigma$-coordinate in which the layers are uniformly distributed in each column to fit the topography. Thus $h$ has no vertical structure in panel c.

For the other panels it is important to remember that the k-index increase downward in the model; k=1 (fortran convention) or k=0 (python convention) is at the surface. Panel d has a region of uniform resolution (~100m) at low k which transitions to vanished thickness (~0) at some value of k in each column. You can sort of see the topography (blue region) upside down. Panels a and b look similar and have a lot of vanished regions with even surface layers being vanished.

Before making sense of these thickness distributions, let's check that the total thickness in each column looks like the topography:

We see that although the thickness distributions are quite different between each model, the total thickness of each column reflects the topography we plotted earlier based on the ocean_geometry.nc file.

The layer thickness almost provides enough information to calculate the actual position of quantities with which we could then make plots. The missing information is the absolute position of the top or bottom.

Interfaces delineate, or bound, layers. The thickness of a layer, $h_{i,j,k}$, is related to the absolute position of the interface above, $z_{i,j,k-\frac{1}{2}}$, and below, $z_{i,j,k+\frac{1}{2}}$ by $$ h_{i,j,k} = z_{i,j,k-\frac{1}{2}} - z_{i,j,k+\frac{1}{2}} \;\;\; \forall \; k=1,2,\ldots,nk $$

where, by convention, integer-valued indices indicate layer-centered quantities and half-valued indices indicate interface-located quantities. Interface and and layer quantities are thus staggered in the vertical.

The diagnostic variable e is the absolute vertical position of model interfaces. Because half-integer indices are not meaningful in most copmuter languages, there is a offset convention as follows.

In FORTRAN: $$ h(i,j,k) = e(i,j,k) - e(i,j,k+1) \;\;\; \forall \; k=1,2,\ldots,nk $$

where arrays indices normally start at 1.

In python $$ h[k,j,i] = e[k,j,i] - e[k+1,j,i] \;\;\; \forall \; k=0,1,\ldots,nk-1 $$

where array indices start at 0. We have also indicated the [k,j,i] order of indices that arises from reading data from an model-generated netcdf file.

Let's look at where the interfaces are by plotting a line for each interface (note the use of the transpose .T operator to get these lines plotted in the right direction):

You can now begin to discern the nature of the coordinates in each mode. If we zoom in on the shelf-break region we will be able to more clearly see what each coordinate mode is doing.

The $\sigma$-coordinate (c) always has nk layers with finite thickness (uniformly distributed) in each column. The $z$*-coordinate model (d) seems to have a variable number of layers but in fact the layers thicknesses vanish wherever the layer would be below the topography. The isopycnal coordinates, both in layer-mode (a) and ALE-mode (b), have on one thick layer on the shelf and fewer finite-thickness layers off-shelf than the other models. In these cases, there are vanished layers at both the top and bottom of the column.

So now we know the location of the interfaces we can presume the center of the layer is in between at $(e[k,j,i]+e[k+1,j,i])/2$. Let's use contourf to shade salinity at the layer centers. Note how we have to create a 2D "x" coordinate to pass to contourf since contourf expects both coordinate arrays to be 2D if either one of them is 2D. We do this by using an expression x=xh+0*z which uses numpy's "broadcasting" feature (see http://docs.scipy.org/doc/numpy-1.10.1/user/basics.broadcasting.html for explanation of rules). We will also plot the interface positions on top of the shaded contours:

The above looks closer to what one imagines things look like but there are some very big problems with the above plots.

1) The apparent topography (white regions at bottom of plots) is quite different between the panels. This happens because contourf only shades between cell centers and so only half of the edge cells are plotted. contourf does not extrapolate beyond the coordinate provided for the data location. We loose a half cell of shading at the top and bottom of the column and also at the left and right of the plot. In the isopycnal-like coordinates, the bottom layer is thick and so we loose a lot.

2) The shading within layers and between columns is interpolated which is introducing interior features and gradients which should not be there. The overlaid interface positions make this apparent for the layer mode for which salinity is absolutely constant along a layer (recall first plot of salinity).

To get the plot we want we essentially need to insert a layer of extra data at the top and bottom of the model and an extra column at both ends. To illustrate lets see how one might do this first for the "layer" output. We'll define a little function to help:

So not the data appears to be plotted from the surface down to the topography. Using this approach for all the coordinates:

A remaining issue is why does there appear to be a salinity inversion in the z*-coordinate model. Technically there is (see salinity plots in i,k-space) but the layers are vanished so we should not be seeing them. This is because contourf is interpolating between thick and vanished layers. The bottom-line is that contourf is assuming the data is smooth and interpreting data inconsistent with the model formulation which considers the data to be piecewise.

Use pcolormesh() to visualize

The most consistent tool for visualizing piecewise data is pcolormesh. An important distinction between contourf and pcolormesh is that the latter takes the coordinates of the corners of cells when shading cell-centered values. Recall we loaded the coordinate yq and inserted an extra value on the left edge - we'll use that for the horizontal coordinate of cell edges. We will horizontally average to get an approximate position for the cell corner heights:

In the above plots, the vanished layers are reasonably hidden and the overall shading for salinity more similar between the plots.

The above method treats each cell as a trapezoid with corners shared between neighboring cells. It does not preserve the mean depth of the cell boundaries. To give more faithful rendering of the what the model can do, a tool is provided in MOM6-examples/tools/analysis/m6toolbox that returns arguments that can be passed straight to pcolormesh consistent with various interpretations of the grid structure, e.g. pcm (piecewise constant thicknesses), plm (piecewise linear), linear (as described above).